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Abstract. Granular materials present memory effects when submitted to tapping 
processes. These effects have been observed experimentally and are discussed here 
in the context of a general kind of model systems for compaction formulated at a 
mesoscopic level. The theoretical predictions qualitatively agree with the experimental 
results. As an example, a particular simple model is used for detailed calculations. 
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1. Introduction 

A granular material is a system composed of a "large" number of macroscopic particles 
or grains whose interactions are inelastic, so that mechanical energy is not conserved. 
Here, in practice, large means in many cases a few thousands or even hundreds of 
particles, as opposite to usual molecular systems having a number of particles of the 
order of Avogadro's number. The macroscopic size of the particles implies that the 
usual concepts of thermodynamics can not be directly translated to these systems. For 
instance, if we consider a typical granular system such as sand, the energy needed to lift 
a grain by one diameter is more than ten orders of magnitude larger than the thermal 
energy of the grain at room temperature. Moreover, due to the inelasticity of collisions 
it is necessary to supply external (nonthermal) energy to the system in order to generate 
a steady state or to study the relaxation of the system towards a stable configuration. 
Therefore, although the concept of granular temperature is often used in the literature, it 
must be understood just as a measure of the velocity fluctuations in the system, without 
being related with any underlying idea of thermal equilibrium. 

The phenomenology of granular media is very rich, showing many characteristic 
complex behaviours [|I|]. Here we will focus in one them, namely compaction, that can 
be roughly defined as the nonthermal relaxation of a loosely packed system of many 
grains under vertical mechanical tapping or vibration. This problem is of fundamental 
importance to many industrial applications and also raises fundamental theoretical 
questions. Therefore, it is not surprising that in the last years systematic experimental 
investigations and theoretical approaches are being developed to described the dynamics 
of compaction processes, as well as the nature of the final state reached by the system. 

The evolution of the density in vibrated granular materials was investigated in a 
pioneering paper by Knight et al P, ||. The experimental set-up they used consisted 
of monodisperse spherical glass beads in a long thin cylinder mounted vertically on a 
vibration exciter. The packing fraction of the beads was measured by means of four 
capacitors mounted at different heights along the tube (see Figure |l|). The shaking 
intensity T, determined by the maximum applied acceleration normalized by gravity, 
was controlled. The system was prepared in a low density initial state before being 
submitted to a sequence of single shakes or "taps" . Time between taps was long enough 
to allow the system to come to rest, so they were completely independent and internal 
resonances were avoided. It was observed that the density increased monotonically, 
tending eventually to a steady value. The number of taps required to reach the final 
density was very large, often larger than 10 5 . 

The authors found that a good description of the experimental data was obtained 
with an inverse logarithm four parameter fit of the form 

p{t) = p ~- TTbWTIJq ' (1) 

where time is measured in number of taps and the parameters p^, Apoo, B, and t c are 
constants that depend only on the tapping strength T. The above inverse logarithmic 
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Figure 1. Sketch of the experimental set-up used in Reference The density is 
measured at four different heights by means of the capacitors CI — C4. 



expression fitted the data better than other more usual laws such as a single exponential, 
a combination of two exponentials, a power law, or a stretched exponential. Two features 
of the experimental results that are relevant for the later discussion here are 

(i) The slope of the relaxation curve is smaller for smaller vibration intensity T, i.e. 
the relaxation is slower for smaller T. 

(ii) The final steady density is a monotonic decreasing function of T [|J. 

A theoretical attempt to formulate a "thermodynamic" description of the steady states 
reached by the system in the long time limit has been carried out by Edwards and 
coworkers @, • They extended the methods and concepts of usual statistical physics 
to granular media. The basic idea is that the volume plays in these systems a role 
analogous to that played by the energy in molecular (elastic) systems. Although there is 
no experimental verification of this theory up to now, it has been shown to be consistent 
with the behaviour of some simple models for granular compaction || [|, [Tcfl . 

In addition to the slow relaxation described above, vibrated granular materials 
exhibit a series of properties that are reminiscent of the typical behaviour of conventional 
structural glasses. This includes effects as annealing, i.e. slow "cooling" properties, and 
hysteresis, when submitted to processes in which the tapping intensity is monotonically 
increased and decreased 0, [Ll|] . This apparent glassy nature of granular compaction led 
Josserand et al to investigate the response of a vibrated granular system to sudden 
perturbations of the intensity V. Their work was suggested by classical experiments 
for the study of aging in glasses, and the realization that the vibration intensity plays 
in vibrated granular media a role analogous to the temperature of molecular systems. 
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Previously, a similar process had been studied by means of numerical simulations in a 



model for compaction |T3|| . 



In the simplest compaction experiment of Reference [12], the vibration intensity was 
instantaneously changed from a value Ti to another T 2 after t w taps. For r x > T 2 it was 
observed that on short times scales the compaction rate increases, while for r 2 > I\ the 
system dilates for short times. Both results are opposite from the behaviour at constant 
T, as discussed above. After several taps, the "normal" behaviour was recovered. This 
is a direct evidence of the presence of short-term memory effects in the system, so that 
the density for t > t w is not determined by its value at t = t w . Other experiments are 
also reported by the authors showing the same kind of non-Markovian behaviour. In 
the next section we will investigate whether this anomalous response can be understood 
in terms of simple and general arguments. 

2. Mesoscopic description of the density evolution 

When a granular system is being vibrated, different kind of events take place in the 
system. The existence of a steady density for constant vibration intensity suggests the 
presence of two different kinds of elementary processes, the ones trying to increase the 
density and the others trying to decrease it. Stationarity arises when both tendencies 
cancel each other. Then, we will assume that the time evolution of the density p in 
a tapping process can be described at a mesoscopic level by a balance equation of the 
form 

^ = /i(r)^)-/ 2 (i> 2 (t), (2) 

where t is measured in number of taps, and the several functions verify 

A(r)>o, / 2 (r)>o, (3) 
mit) > o, /i 2 (t) > o. (4) 

The quantities pi(t) and p 2 {t) are assumed to include all the correlation effects. Note 
that we do not assume that p\ and p 2 can be expressed as functions of p(t), so that 
Equation (Q) is not closed. Our aim in the following will be to investigate what 
information can be obtained about fi, f 2 , p±, and p 2 by using very general and plausible 
physical arguments. First, we notice that since we are measuring the time in units of 
complete taps, if there is no tapping there is no time evolution either. Therefore, it 
must be 

A(r = o) = / 2 (r = o) = o. (5) 

Moreover, the number of elementary processes, both tending to increase and to decrease 
the density, are expected to increase with T. This expectation is reinforced by the 
analogy of T with the temperature, as mentioned in the previous section. Then we 
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assume that /i(r) and f 2 (X) are both monotonic increasing functions. Let us introduce 
the function 

rfr)-£g. (6) 

representing the ratio of the rate associated to the decompaction processes to that of the 
compaction ones. As already mentioned, experiments show that the steady density is a 
decreasing function of T. This suggests that g(T) must be an increasing function of the 
vibration intensity, so that decompaction events become relatively more relevant. This 
assumption can be further justified by analysing the behaviour of the steady solutions 
of Equation (Q) [Tl|. Then we can rewrite our mesoscopic evolution equation as 

^p(t) = f 1 (T)[p 1 (t)-g(T)p 2 (t)}. (7) 

Here /i(r), g(T), pi(t), and p, 2 {t) are an positive quantities. Moreover, the first two 
ones are increasing functions of T, vanishing in the limit r — > 0. This fully specifies 
the general kind of evolution equations we will deal with. Let us consider the following 
experiment. A system is vibrated with a constant intensity T. At a given time t w , the 
intensity is instantaneously changed to T + Ar. We want to study the change in the 
relaxation rate r(t) = dp(t)/dt. Just before the change we get from (0), 

r w = f 1 (T)[ f i 1 (t-)-g(T)fi 2 (t-)], (8) 
while just after the change it is 

r' w = h{V + Ar) [ Ml (t+) - g(T + Al> 2 (t+)] . (9) 
Then, for Ar < 1, we get 

^ = ! ^? £ = (10) 

where we have defined the function (r(t) and p. 2 {t) are computed over the relaxation 
curve at constant T) 

m = J^r(t)-h(T)g'(T)p 2 (t). (11) 

Upon deriving Eq. (|T0| ) we have assumed that pi(t) and p 2 (t) are continuous at t = t w . 
This seems physically plausible since they are functionals of the state of the system. 
For instance, in the mesoscopic description they can be expressed as given moments of 
the underlying distribution function of the system, which is known to be a continuous 
function of time. In the limit t w — > oo, the system has time to reach the stationary 
density corresponding to the intensity T before the change, so that r w = and \(t w ) < 0, 
because of the properties of the functions /i(r) and g(T). In the opposite limit, t w — > 0, 
and assuming that the initial density is close to its minimum value, it is X(t w ) > 0. To 
derive this, we have taken into account that p 2 must vanish in the lowest density limit, 
where by definition no processes decreasing the density are possible. The conclusion of 
this discussion is that \(t w ) has opposite signs in the short and large t w limits. Then, it 
follows from Eq. ([ToD that when the system has been vibrated for a short time before the 
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change of vibration intensity, Ar w and Ar have the same sign, i.e. the response of the 
system is what we have described in the previous section as normal. On the other hand, 
if the vibration time before the change is large, the response of the system is anomalous. 
In this context, the experiments reported in Reference [E| would correspond to "large" 



vibrating periods before the abrupt change in the shaking intensity. 
3. A simple model for granular compaction 

In this Section we are going to describe a one-dimensional model for granular compaction 
[15| that is simple enough as to allow for detailed calculations, and at the same time 



captures many of the experimental characteristic features. We consider a lattice in which 
each site i can be either occupied by a particle or empty (occupied by a hole) . A variable 
rrti is assigned to each site, taking the value = 1 if the site is empty, and rrii — if 
there is a particle on it. The time evolution of the system is defined in the following 
way. The only possible elementary events occurring in the system are the adsorption of 
a particle on an empty site and the desorption of a particle from the lattice to the bulk. 
Then the dynamics is formulated by means of a master equation with a transition rate 
for the change of rrii into 1 — given by 

Wi(m) = ^(m 4 _! + m i+1 ) [e + rm{\ - 2e)\ . (12) 

Here m = {mi,m2, . . .}, v is a frequency defining the characteristic time scale, and e is 
a dimensionless parameter taking values in the interval < e < 1. Thus the transition 
rate for the adsorption of a particle at site i is 

W+(m) = V -{\ -z)mi{mi- x +m i+1 ), (13) 

and that for the desorption of a particle 

W-(m) = ^e(l - mi){mi-i + m i+1 ). (14) 

It is seen that for e = 1 no particle can be absorbed, while for e = desorption processes 
are inhibited. Moreover, all processes are restricted, in the sense that a site can change 
its state only if at least one of its nearest neighbour sites is empty. A similar kind of 
facilitated dynamics has been considered in the formulation of Ising models for glassy 
relaxation [|16|]. A physical picture of the model can be obtained by associating a hole 
with a region of the granular system having a low density, and a site occupied by a 
particle with a high density region. The facilitated dynamics tries to translate the idea 
that fluctuations leading to rearrangements in a region require a neighbour region with 
low packing fraction. Our one-dimensional model is a very idealized representation of a 
low horizontal layer in a vibrated granular system, the parameter e being related with 
the intensity of vibration. 

In order to model a tapping experiment, the system is initially placed in a purely 
random configuration, from which it relaxes with e = until it gets eventually trapped 
in a metastable state, with all the holes surrounded by two particles. This is a low 
density configuration that is taken to correspond to the loosely initial conditions used 
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Figure 2. Time evolution of the density of the model described in the main text when 
it is submitted to a tapping process with e = 0.5. Data obtained with three different 
values of to, namely 0.002, 0.001, and 0.02 have been plotted. 



in real experiments. Then, a pulse is generated by instantaneously increasing e to a 
given value for a time period t^. Afterwards, the system relaxes again with no external 
excitation, i.e. with e = 0. In order to mimic what is done in experiments, the relaxation 
lasts long enough as to allow the system to reach a metastable configuration, from which 
no evolution is possible in absence of external excitation. This completes a tap event. 
The process is repeated to generate a series of them. 

A first test of the relevance of the model to study compaction in vibrated granular 
systems is, of course, whether it leads to the same kind of density evolution at constant 
vibration intensity as experimentally observed. We have verified that this is the case 
in the limit euto <C 1. The mean density in an homogenous configuration is given by 
p(t) = 1 — (mi(t)), with the angular brackets denoting ensemble average and i being 
arbitrary. 

An example is given in Figure |^, where the relaxation of the density is plotted as 
a function of the scaled time r n = vet^n, n being the number of taps. The numerical 
data have been obtained from Monte Carlo simulations and different values of e and 
to have been used. The fact that all curves collapse indicates that r n is the relevant 
time scale for the compaction problem. Also plotted in the same figure is the fit to 
the phenomenological law (]]]), with parameters p^ = 1.10, Apoo = 0.40, B = 0.39, 
and t c = 3.37/^eto- It is seen that the inverse logarithmic law describes very well the 
simulation results. A similar behaviour has also been found in other models for granular 



compaction [17|, 0, [L9| . Nevertheless, we ought to say that we have not been able to 
derive the heuristic relaxation law by analytical methods, in spite of the tractability of 
our model. It is possible that it be just a convenient fitting expression over a wide time 
window. A strong indication supporting this idea is that the steady density predicted 
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Figure 3. The effective dynamics tries to define the transitions between the crosses, 
which correspond to the results of the successive taps. 



by the logarithmic law in the limit of an asymptotically large number of taps not 
only is in disagreement with the simulations, but is clearly unphysical since it is larger 
than one. The same happens with the experimental results, where poo is sometimes 
greater than the random close packing 0. In fact, it is possible to derive an analytical 
expression for the asymptotic density reached by the model with the result 

p M ~ i _ l eut0i (15) 
valid again in the limit euto <C 1. Also this expression has been checked by Monte Carlo 



simulations 15 



4. Effective dynamics for tapping processes 

The stochastic model formulated in the previous section, with the transition rates given 
by Eq. (0), can be used to study a variety of compaction processes by specifying the time 
dependence of the control parameter e. In particular, we have already discussed some 
applications to tapping processes, based on the particularization of the general equations 
following from the master equation to a specific way of vibrating the system. A different, 
and more appealing, approach is to look for an effective master equation following from 



(|T2|) and being appropriated for a given experiment. We have developed such a program 
for tapping processes ||. The idea is to look for effective transitions rates, W e /(m|m'), 
connecting the initial and final states of the system when it is submitted to a "elementary 
event" (see Figure [$]). The latter is defined as the combination of a tap and the posterior 
free relaxation to a metastable configuration. In the limit of ez/t <C 1, three groups of 
possible transitions are identified: 

a) Elementary diffusive events, conserving the number of particles. They correspond to 
the interchange of a hole and a particle, 

. ..100... -»• ...010..., ...001... -»• ...010.... (16) 

The effective transition rate for each of these processes is a/2, where 

a = — — . (17) 
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is a positive constant playing a role similar to that played by T in the real experiments. 

b) Transitions increasing the number of particles. There are three of them, 

. ..101... -»• ...010..., (18) 
with transition rate a/2, 

. ..101... -> ...001..., (19) 
with transition rate a/4, and 

...101... -»• ...100..., (20) 
also with transition rate a/4. 

c) Transitions decreasing the number of particles. These are 

...00100... -»• ...01010..., (21) 
with transition rate a 2 /2 and 

...01000... -»• ...01010..., (22) 
...00010... -»• ...01010..., (23) 

both with W e f = a 2 /4. Only those variables corresponding to sites whose state is 
changing or conditioning the transition are represented in the above expressions. The 
equivalence of this model to the original one defined by the transition rates (|i~2|) in 
the limit of small a, has been tested by comparing the Monte Carlo simulation results 
obtained with both models. In fact, the data show that, for the density relaxation, 
the results from the effective model and the original one differ by less than 2 per cent 
for a < 0.5. For constant a ^ 0, the system evolves from the initial low density 
configuration to a final state characterized by a density 

Ps {a)= l -[l + (l + Aa)- l/2 \. (24) 

A detailed analysis of the properties of the steady state || shows that they are 
consistent with Edwards's theory |5j, ||, [7|], with the compactivity being identified as 
X = —(In a)" 1 . Furthermore, when the system described by the effective transition 
rates is submitted to processes in which the tapping intensity is first monotonically 
increased and afterwards decreased also monotonically, its time evolution presents the 
reversible-irreversible branches observed in experiments [[2(J. In Figure f| an example 



of the response of the system to one of these cycles is presented. The compactivity X 
of the system is decreased and increased with the same rate, r = 10~ 5 . Starting from 
the loosest configuration, for large enough vibration intensity (or compactivity) the 
density of the system approaches the steady curve. Afterwards, when the compactivity 
is decreased and increased, always with the same rate, two (approximately) reversible 
curves are obtained. These hysteresis effects are related with the existence of a "normal 
evolution curve", fully determined by the program of increase of the intensity, and 
having the strong property of attracting any other solution of the master equation with 
the same program [ p0|] . 
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Figure 4. Evolution of the density when the system is submitted to a tapping cycle 
as described in the text. The diamonds and the crosses represent the approximately 
reversible cooling and heating processes, respectively. 



5. Memory effects 

Now, we will show that the effective dynamics introduced above leads to a model 
for tapping processes that belongs to the general class discussed in Section || and, 
consequently, also presents the short term memory effects observed in vibrated granular 
materials. We start by noting that the steady density for constant intensity, given by 
is a monotonic decreasing function of a, going from p max = 1 to p m j„ = 0.5. Next, 



from the master equation with the transition rates (|Tq ) -(|23 D it is obtained 



dp a 

— = axiQi(t) 

dt K J 2 



2 r 



£ooioo(£) + -^oiooo (t) + - x oooio(t) 



(25) 



where xxoi is the concentration of hole-particle-hole clusters, xooioo is the concentration 
of two particles-hole-two particles clusters, and so on. Comparison of the above equation 
with (0) shows that both equations have the same form, with the choices 

fi(a) = g(a) = a, (26) 

in(t) = xioi(t), nz(t) = ^xooiooOO + \ koioooOO + £00010 (*)] • (27) 

Since f\ and g are monotonic increasing functions of ck, and /^i and /i 2 are defined positive 
quantities, the model verifies the conditions required by the validity of the discussion 
carried out in Section |^. Therefore, we can directly write from Eq. fllOD 

^ = A(U (28) 

where the function X(t) determining the nature of the response of the system is 

r(t) 

X(t) = -ap 2 {t). (29) 
a 
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Figure 5. Response of the system to a sudden change in the vibration intensity from 
a = 0.15 to a = 0.125. When the jump is introduced at t w = 10 < t c , the response 
is normal, i.e. the slope of the relaxation curve decreases. On the other hand, for 
t w — 50 > t c an anomalous response is observed. 



For instance, for a = 0.15, Monte Carlo simulations show that \(t) > for t < t c , while 
X(t) < for t > t c , with 19 < t c < 20 |14|]. According to the theory presented here, 



when the vibration intensity a is modified at a time t w < t c , a normal response in which 
the intensity jump and the relaxation rate jump have the same sign is to be expected. 
On the other hand, when t w > t c , stimulus and response should have opposite signs, 
corresponding to an anomalous response. In Figure [5] we have plotted the time evolution 
of the density of a system which is being vibrated with a = 0.15, and the intensity is 
suddenly decreased to a = 0.125 at t — t w . The only difference between the two curves 
is in the value of t w ; in one case it is t w = 10 < t c while in the other t w = 50 > t c . 
In agreement with the theoretical analysis, the relaxation rate at the jump decreases in 
the first case and increases in the latter. 

It is important to stress the generality of the arguments in Section ^. Although we 
have restricted ourselves in this work to a particular simple model of compaction, the 
theoretical scheme presented there is rather general. For instance, an equation like (0) 



is also found for the so-called parking model [19], ^T[|. Even more, the one- dimensional 
Ising model with Glauber dynamics also belongs to the group | 22fl . In summary, the 
memory effects discussed in the context of compaction in granular materials seem to be 
quite general effects showing up in many different systems. 
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